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Abstract 

In this paper a boundary integral numerical method for the dynamics of non-spherical cav- 
itation bubbles in inviscid incompressible liquids is described. A very efficient feature of 
this approach is that it involves only surface values of the velocity potential and its first 
derivatives, thus avoiding the problem of solving the Laplace equation in the entire domain 
occupied by the liquid. To i3 lustrate the performance of the method the collapse of a bubble 
in the vicinity of a solid wall and the collapse of three bubbles with cc: linear centers are 
considered* 


Introduction 

The development of the study of the dynamics of non-spherical bubbles has been seriously 
hindered in the past by the lack of computational tools applicable to severe deformations of 
the spherica'*. shape. Existing methods **** are cumbersome to implement and require long exe- 
cution times. These problems have prevented their widespread usage and common acceptance. In 
this study we shall present a much more efficient method based on the boundary integ^^al ap- 
proach ^ which can be implemented in a relatively small number of program instructions (500- 
600). It is flexible, precise, and requires short execution times (typically 3-5 minutes on 
a Univac 1100). The method is applicable to irrotational , incompressible, inviscid flow and 
we present an axisymmetric version of it. To illustrate the results obtainable in this way 
we consider the Rayleigh problem, one of the cases studied by Plesset and Chapman ^ , and the 
collapse of three bubbles with collinear centers. 

Mathematical formulation 


We consider the flow induced by one or more bubbles in a liquid occupying a domain 
bounded by the surface of the bubbles, S , rigid boundaries, S , and a "surface at infinity" 
S^. In the hypotheses of inviscid irrotational flow the velocity field u can be derived from 
a potential, u = V<t>, and the mathematical problem can be put in the form^ 


= 0 in n , 

= 0 on S , 

5n r 

d(J) 1 , i j ^b 

5 ? on 

dx 

— - 7$ on , 

(j, - 0 , p - on , 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


where p is the pressure, p is the density, 3/3n is the normal derivative, and d/dt=3/3t+V^» 
is the convective der ivative.The pressure in the bubble, p^, and at infinity, p , are taken 
constant in this study, although this is by no means required by the method* Surface tension 
effects have been disregarded, but they can be included with relative ease* 


The strategy of solution is in principle quite straightforward. Suppose that the position 
of the bubble surfaces Sj^ and the value of ^ on Sj^ are kno%m at time t* Then Eq*(1) can be 
solved using the boundary conditions (2), (5), and the known ^ on From a knowledge of the 
potential it is then possible to evaluate the right-hand sides of (3) and (4), which allows 
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the computation of the new configuration of and of the potential on it at time t+At. The 
cycle is then repeated. 

The crucial step in this procedure is the solution of the potential problem* In the past 
this has been done by finite differences or by singularity methods The first ap- 

proach is wasteful because it requires that: the potential be known over the entire domain U, 
although only free surface values are needed . In addition t the position of does not co- 
incide with grid nodes in general , which leads to well-known problems of accuracy » coding, 
and others. The second method avoids the explicit solution of (1) because the assumed form 
for 0 already is a solution of the Laplace equation. However a certain amount of guessing is 
necessary, and it is not clear whether very extreme bubble deformations can be described 
satisfactorily with this approach. Our technique, based on the boundary integral, or bound- 
ary element, method^, retains the advantages of the singularity method while avoiding its dis- 
advantages since it is based on an exact relation. Green’s Identity. This identity for a 
point X belonging to or has the form* 

° jnx-x’l~* If,- ■l'(x*) I;;. |x- 2 E'l~‘) ds* . (6) 

Sb+Sr 

In conditions of axial symmetry, which we assume, this relation can be simplified by carrying 
out the integration over the azimuthal angle explicitly with the result 
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In this equation x and x' belong to a meridian (half) plane passing through the axis of sym- 
metry and Tjj and are the traces on this plane of and respectively; s' is the arc 
length on these curves F. In the meridian plane we choose a Cartesian system of coordinates 
(r,z), of which the z-axis coincides with the axis of symmetry- Then G,H have the form 

G(r,z;r',z') = K(m) , H(r,z;r'z') K(m) /A} , (8) 

7t/A tt an 

where A ~ (r+r')"+ (z-z')^, m=4rr'/A, and K denotes the complete elliptic integral of the 
first kind. 


In the applications to be described in this paper the existence of the rigid boundary is 
accounted for by introducing a system of "image" bubbles. Therefore condition (2) will be 
implicitly satisfied and it is no longer necessary to indicate explicitly the rigid boundary 
Sj. or Fj.. Accordingly, we shall not carry these symbols along in the following equation^ and 
we shall denote simply by F . 


Numerical method 


Since 4> is known on F from the previous time step, Eq.(7) can be regarded as an integral 
equation for ?(t/3n'. To solve it, we approximate the line F by a polygonal line consisting 
of n segments Fj on each of which and34>/3n' are taken constants. (We have also made some 
numerical experiments assuming a linear variation of obtaining virtually identical results) . 
In this way we obtain from (7) a linear system 
n 
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where w * 34>(y.)/ 3n' can be 
point yj of Fj ^and 
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interpreted as the value of the normal velocity at the mid- 
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Here to be interpreted as the coordinates of the midpoint ^ of r . The integrals 

in (10) are performed numerically by a six-point Gaussian formula except for those over 
for which the integrand has an (i*^^egrable) logarithmic singularity. A series expansion in 
the neighborhood of this singularity coupled with Simpson's formula over the rest of the 
interval is used in this case. The elliptic integrals are computed using the simple approxim- 
ate formulae of Ref. 8. The tangential velocity along the bubble surface Is computed to the 
same accuracy as the normal velocity as follows 


v(x.) 


iilTrilT 


( 11 ) 


where x^+^ and x^^ are the extrema of the segment The normal and tangential velocities 
computed in this way are then referred to the (r,z) system of coordinates and transported 
by linear Interpolation from the midpoints to the extrema of the segments r^. 

At this point Eqs.O) and (4) can be used to compute the advanced-time values of the 
positions of the vertices of the polygonal line and of the associated velocity potentials. 
For this purpose the following second-order formula is used 
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where x» At / At and At =t ^^-t .A similar relation is used to compute the new values 
*^i' ^4 positions of tRc possibility of using higher order time integration formu- 

lae of this type is a distinctive advantage of the present method over finite-differences 
ones. This advantage stems from the fact that here it is possible to follow the individual 
trajectories o.‘ ♦^he vertices of the polygonal line during the entire calculation. The time 
step used in (12) is adjusted during the calculation so as to prevent excessive variations 
of 2 ^ during a single step. 


Results 




We shall present the results in terms of lengths and times made dimensionless with re- 
spect to Rq, the initial bubble radius^ and t^-R^ { p/ (Pa,“Pb) ^ • The pressure inside the 
bubble and at infinity have been kept constant in all the examples discussed. 

To test the reliability of the method and of the code we have first of all computed the 
collapse of a single spherical bubble in an unbounded liquid, the well-known Rayleigh prob- 
lem . The calculation started to develop instabilities in the velocity for a bubble radius 
R==1.411 X 10"“^, radial velocity R=483.7, at time t=0.9162. It is well known that this spher- 
ical collapse is unstable . Therefore, the ability of the code to compute a decrease in rad- 
ius by nearly two orders of magnitude with no smoothing techniques applied and within a max- 
imum 1.8% deviation from sphericity is an indication of an excellent performance. The ana- 
lytical result for the radial velocity for the value of the radius indicated above is 
R«487.4, again in excellent agreement with the numerical result. Finally, the analytical 
value of the total collapse time is tc=0. 91468, in very good agreement with the computed 
value, which practically corresponds to total collapse. 

As a second test we have considered one of the cases studied by Plesset and Chapman^ .Here 
the bubble collapses in the presence of a plane rigid wall, from which it is separated by a 
distance equal to half the initial radius at t>0. In the calculation the plane was simulated 
by introducing an image bubble equal and symmetrically located with respect to the real bub- 
ble. By taking advantage of the symmetry of the problem it is possible to <aaintain the num- 
ber of unknowns in the system (9) equal to the number of segments of thcs real bubble. We 
show in Fig.1 the bubble shapes at selected instants (heavy lines) and the trajectories of 
some of the points on the bubble surface (light lines). It is interesting to observe the very 
strong convergence of streamlines in the jet which evolves in the l«iter stages of the collapse. 
The shapes of Fig.1 agree with those of Ref.1, but a more stringent comparison is afforded by 
the values of the velocity of the "north pole" of the bubble as computed by us (Fig. 2, upper 
line) and as given in Ref *1 (Fig. 2, open circles) . Although the comparison is good, some minor 
discrepancies exist the origin of which is not clear at the present time. It may be noted how- 
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ever from the shapes published by Plesset and Chapman ‘ that some imprecision affects their 
method near the ** north pole** of the bubble^ where the tangent falls to be zero as It should. 
The lower line of Fig. 2 shows the time development of the velocity of the "south pole" of 
the bubble. 

As a final example we consider the collapse of a system of three equal and equally spaced 
bubbles, with centers on the axis of symmetry. The evolution of the process is shc%m in Pig. 

3 for initial bubble spacings of 0.5 (left), 0.2 (center), and 0.1 (right). Because of the 
symmetry of this situation only the upper bubble and the upper half of the middle bubble are 
shown in the figure. The collapse of the central bubble is strongly inhibited by that of the 
other two and its "poles" move very little. This effect of course increases with the proxim- 
ity oi the bubbles. 

Finally, to give an idea of the computational requirements we present in Table 1 some in- 
formation on the number of segments used (for all the examples this is the number of segments 
in the first quadrant of the (r,z) plane) . the number of time steps, and the CPU needed for 
the calculation on a UNIVAC 1100/8. Note that we have not tried yet to optimize the method 
with respect to time step size and number of segments used. For instance, it is believed that 
the latter quantity, in the cases of Fig, 3, could be reduced substantially without appreci- 
able loss of accuracy. 


Table 1 . Some details on the computation 


Case 

Number of segments 

Number of time steps 

CPU time (mins and secs) 

Rayleigh problem 

16 

105 


7 

Figure 1 

32 

79 

3 

15 

Figure 3, left 

48 

92 

8 

30 

Figure 3, center 

48 

99 

9 

12 

Figure 3, right 

48 

100 

9 

7 


Conclusions 


The boundary integral method described in the present study has proven to be a very effi- 
cient and accurate tool for the study of non-spherical bubble dynamics. It makes possible 
extensive investigations of non-spherical bubble behavior at a fraction of the effort and 
cost required by other methods. In addition to its value for specific problems, we think that 
this feature is particularly important in the present state of research In this field because 
in our opinion only numerical experiments can help develop the intuition necessary for fur- 
ther progress. 
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Figure :• Coliapse of an initially spherical bubble in tie neighborhood of a solid plane 
fro' which it is separated by 0.5 times the initial radius. The shapes shown (heavy lines) 
’^^e at t«0, t*0.950, t«1.010, and t«1.C33. The collapse was completed during the time step 
after the last one shown. Light lines are particle trajectories. 






Figure 2. Comparison between the velocity of the "north pole 
of the bubble of Fig, 1 (upper line) and the values given by 
Plesset and Chapman * (open circles). The lower line shows the 
velocity of the "south pole**. 



Figure 3. Collapse of three equal and equally spaced bubbles with coll inear centers. The process 
is symmetric about the z«0 plane. The initial bubble spacings are 0.5 (left) ^ 0.2 (center), and 0. 
(right) tiroes the initial radius. The shapes shown are for: A, t=0; B, t»0.950; C, t=1.063; D, 
t»1.105 (left and right), t=1.103 (center); E, t»1.130; F, t~1.135. 


